#include "coupled_pairwise.hpp"
#include "parcsr_matop.hpp"
#include "MG.hpp"

#define NORMAL
// #define THRESHOLD
// #define TOPK

template<typename idx_t, typename setup_t, int dof>
idx_t * coupled_pairwise(const par_CSRMatrix<idx_t, setup_t, setup_t, dof> * par_fine_mat,
    par_CSRMatrix<idx_t, setup_t, setup_t, 1> * & R_mat, par_CSRMatrix<idx_t, setup_t, setup_t, 1> * & P_mat,
    const setup_t theta)
{
    const MPI_Comm comm = par_fine_mat->comm;
    int my_pid; MPI_Comm_rank(comm, & my_pid);
    int num_procs; MPI_Comm_size(comm, & num_procs);
    constexpr MPI_Datatype mpi_idx_type = sizeof(idx_t) == 8 ? MPI_LONG_LONG : MPI_INT;

    const idx_t my_num_F = par_fine_mat->end_row - par_fine_mat->beg_row;

    par_CSRMatrix<idx_t, setup_t, setup_t, 1> S_mat(comm, par_fine_mat->rows_partition, par_fine_mat->cols_partition);
    S_mat.commpkg = par_fine_mat->commpkg;
    S_mat.col_map_offd = par_fine_mat->col_map_offd;
    assert(S_mat.beg_row == S_mat.beg_col && S_mat.end_row == S_mat.end_col);
    const idx_t _fbeg = S_mat.beg_row, _fend = S_mat.end_row;

    constexpr int e_size = dof*dof;
    auto calc_blk_norm = [] (const setup_t * val) {
        setup_t ret = 0.0;
        for (int f = 0; f < e_size; f++) ret += val[f] * val[f];
        // ret = val[4]*val[4] + val[5]*val[5] + val[7]*val[7] + val[8]*val[8]; assert(dof == 3);
        return sqrt(ret);
    };
    std::vector<idx_t> bdr_F;// 细网格上位于进程边界的点，不含孤立井
    bool sorted = false;
    auto sort_and_find_bdr = [&S_mat, my_num_F, _fbeg, & bdr_F, & sorted] () {
        for (idx_t i = 0; i < my_num_F; i++) {
#define PRES_COL_ORDER
#ifdef PRES_COL_ORDER // 保证列号升序（以便coupled和naive在单进程时严格复现）
            sort_vals_with_indices<idx_t, setup_t>(S_mat.diag.row_ptr[i + 1] - S_mat.diag.row_ptr[i],
                    S_mat.diag.vals    + S_mat.diag.row_ptr[i],
                    S_mat.diag.col_idx + S_mat.diag.row_ptr[i] );
#else // 快排不稳定排序
            my_qsort1<idx_t, setup_t, idx_t, -1>(   S_mat.diag.vals    + S_mat.diag.row_ptr[i],
                                                    S_mat.diag.col_idx + S_mat.diag.row_ptr[i],
                                        0, S_mat.diag.row_ptr[i + 1] - S_mat.diag.row_ptr[i] - 1);
#endif
            if (S_mat.offd.row_ptr[i + 1] > S_mat.offd.row_ptr[i]) {
#ifdef PRES_COL_ORDER // 保证列号升序（以便coupled和naive在单进程时严格复现）
                sort_vals_with_indices<idx_t, setup_t>(S_mat.offd.row_ptr[i + 1] - S_mat.offd.row_ptr[i],
                    S_mat.offd.vals    + S_mat.offd.row_ptr[i],
                    S_mat.offd.col_idx + S_mat.offd.row_ptr[i] );
#else // 快排不稳定排序
                my_qsort1<idx_t, setup_t, idx_t, -1>(   S_mat.offd.vals    + S_mat.offd.row_ptr[i],
                                                        S_mat.offd.col_idx + S_mat.offd.row_ptr[i],
                                            0, S_mat.offd.row_ptr[i + 1] - S_mat.offd.row_ptr[i] - 1);
#endif
                // 如果该点是井，则不算入边界点，因为它们不参与聚合，直接作为粗点
                if (iso_wells.find(i + _fbeg) != iso_wells.end()) continue;
                else bdr_F.push_back(i);// 记录这个点为处于进程边界的点
            }
        }
        sorted = true;
    };

    // double t1 = wall_time();
    
    S_mat.diag.nrows = par_fine_mat->diag.nrows; S_mat.offd.nrows = par_fine_mat->offd.nrows;
    S_mat.diag.ncols = par_fine_mat->diag.ncols; S_mat.offd.ncols = par_fine_mat->offd.ncols;
    S_mat.diag.nnz   = par_fine_mat->diag.nnz  ; S_mat.offd.nnz   = par_fine_mat->offd.nnz  ;
    S_mat.diag.alloc_mem(false);                 S_mat.offd.alloc_mem(false);
    for (idx_t i = 0; i <= my_num_F; i++) {
        S_mat.diag.row_ptr[i] = par_fine_mat->diag.row_ptr[i];
        S_mat.offd.row_ptr[i] = par_fine_mat->offd.row_ptr[i];
    }
#if defined(THRESHOLD)
    setup_t * my_dvs = new setup_t [my_num_F];
    for (idx_t i = 0; i < my_num_F; i++) {// 逐行计算Frobenius
        setup_t dv = 0.0;
        for (idx_t p = S_mat.diag.row_ptr[i]; p < S_mat.diag.row_ptr[i + 1]; p ++) {
            S_mat.diag.col_idx[p] = par_fine_mat->diag.col_idx[p];
            S_mat.diag.vals   [p] = calc_blk_norm(par_fine_mat->diag.vals + p * e_size);
            if (S_mat.diag.col_idx[p] == i) dv = S_mat.diag.vals[p];
        }
        for (idx_t p = S_mat.offd.row_ptr[i]; p < S_mat.offd.row_ptr[i + 1]; p ++) {
            S_mat.offd.col_idx[p] = par_fine_mat->offd.col_idx[p];
            S_mat.offd.vals   [p] = calc_blk_norm(par_fine_mat->offd.vals + p * e_size);
        }
        assert(dv != 0.0);
        my_dvs[i] = dv;
        for (idx_t p = S_mat.diag.row_ptr[i]; p < S_mat.diag.row_ptr[i + 1]; p ++) {
            S_mat.diag.vals[p] = (S_mat.diag.vals[p] * S_mat.diag.vals[p]) / dv;// 除以对角线
        }
        for (idx_t p = S_mat.offd.row_ptr[i]; p < S_mat.offd.row_ptr[i + 1]; p ++) {
            S_mat.offd.vals[p] = (S_mat.offd.vals[p] * S_mat.offd.vals[p]) / dv;
        }
    }
    
    idx_t * shrk_diag_rpt = new idx_t [my_num_F + 1]; shrk_diag_rpt[0] = 0;
    idx_t * shrk_offd_rpt = new idx_t [my_num_F + 1]; shrk_offd_rpt[0] = 0;
    const setup_t threshold = theta * theta;
    {
        setup_t * send_buf = new setup_t [S_mat.commpkg->send_map_starts[S_mat.commpkg->num_sends]];
        setup_t * recv_buf = new setup_t [S_mat.commpkg->recv_vec_starts[S_mat.commpkg->num_recvs]];
        constexpr MPI_Datatype etype = sizeof(setup_t) == 8 ? MPI_DOUBLE : MPI_FLOAT;
        for (idx_t s = 0; s < S_mat.commpkg->num_sends; s ++) {
            setup_t * buf = send_buf + S_mat.commpkg->send_map_starts[s];
            const idx_t num = S_mat.commpkg->send_map_starts[s + 1] - S_mat.commpkg->send_map_starts[s];
            const idx_t * loc_ids = S_mat.commpkg->send_map_elmts + S_mat.commpkg->send_map_starts[s];
            for (idx_t k = 0; k < num; k++)
                buf[k] = my_dvs[loc_ids[k]];
            MPI_Isend(buf, num, etype, S_mat.commpkg->send_pids[s], 77, S_mat.commpkg->comm, & S_mat.commpkg->send_reqs[s]);
        }
        for (idx_t r = 0; r < S_mat.commpkg->num_recvs; r ++) {
            setup_t * buf = recv_buf + S_mat.commpkg->recv_vec_starts[r];
            const idx_t num = S_mat.commpkg->recv_vec_starts[r + 1] - S_mat.commpkg->recv_vec_starts[r];
            MPI_Irecv(buf, num, etype, S_mat.commpkg->recv_pids[r], 77, S_mat.commpkg->comm, & S_mat.commpkg->recv_reqs[r]);
        }
        MPI_Waitall(S_mat.commpkg->num_recvs, S_mat.commpkg->recv_reqs, MPI_STATUSES_IGNORE);
        MPI_Waitall(S_mat.commpkg->num_sends, S_mat.commpkg->send_reqs, MPI_STATUSES_IGNORE);

        for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi ++) {
            idx_t dz_cnt = 0, oz_cnt = 0;
            for (idx_t p = S_mat.diag.row_ptr[loc_fi]; p < S_mat.diag.row_ptr[loc_fi + 1]; p ++) {
                const idx_t loc_fj = S_mat.diag.col_idx[p];
                S_mat.diag.vals[p] /= my_dvs[loc_fj];
                if (S_mat.diag.vals[p] < threshold) {
                    S_mat.diag.col_idx[p] = -1;;// 不超过强连通阈值，标记
                    dz_cnt ++;
                }
            }
            for (idx_t p = S_mat.offd.row_ptr[loc_fi]; p < S_mat.offd.row_ptr[loc_fi + 1]; p ++) {
                const idx_t loc_fj = S_mat.offd.col_idx[p];
                S_mat.offd.vals[p] /= recv_buf[loc_fj];
                if (S_mat.offd.vals[p] < threshold) {
                    S_mat.offd.col_idx[p] = -1;// 不超过强连通阈值，标记
                    oz_cnt ++;
                }
            }
            shrk_diag_rpt[loc_fi + 1] = S_mat.diag.row_ptr[loc_fi + 1] - S_mat.diag.row_ptr[loc_fi] - dz_cnt;
            shrk_offd_rpt[loc_fi + 1] = S_mat.offd.row_ptr[loc_fi + 1] - S_mat.offd.row_ptr[loc_fi] - oz_cnt;
        }
        for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi ++) {
            shrk_diag_rpt[loc_fi + 1] += shrk_diag_rpt[loc_fi];
            shrk_offd_rpt[loc_fi + 1] += shrk_offd_rpt[loc_fi];
        }
        delete [] send_buf; send_buf = nullptr;
        delete [] recv_buf; recv_buf = nullptr;
    }
    {
        // long long AS_nnz[2] = {S_mat.diag.nnz + S_mat.offd.nnz, shrk_diag_rpt[my_num_F] + shrk_offd_rpt[my_num_F]};
        // long long GL_nnz[2];
        // MPI_Reduce(AS_nnz, GL_nnz, 2, MPI_LONG_LONG, MPI_SUM, 0, comm);
        // if (my_pid == 0) printf("A %lld => S %lld (%.4f)\n", GL_nnz[0], GL_nnz[1], (double)GL_nnz[1]/(double)GL_nnz[0]);

        S_mat.diag.nnz = shrk_diag_rpt[my_num_F];
        S_mat.offd.nnz = shrk_offd_rpt[my_num_F];
        idx_t * shrk_diag_cid = new idx_t [S_mat.diag.nnz]; setup_t * shrk_diag_val = new setup_t [S_mat.diag.nnz];
        idx_t * shrk_offd_cid = new idx_t [S_mat.offd.nnz]; setup_t * shrk_offd_val = new setup_t [S_mat.offd.nnz];
        std::unordered_map<idx_t, std::pair<int, idx_t> > glb2pidloc;
        for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi ++) {
            idx_t d_cnt = 0, o_cnt = 0;
            for (idx_t p = S_mat.diag.row_ptr[loc_fi]; p < S_mat.diag.row_ptr[loc_fi + 1]; p ++) {
                if (S_mat.diag.col_idx[p] != -1) {
                    shrk_diag_cid[shrk_diag_rpt[loc_fi] + d_cnt] = S_mat.diag.col_idx[p];
                    shrk_diag_val[shrk_diag_rpt[loc_fi] + d_cnt] = S_mat.diag.vals   [p];
                    d_cnt ++;
                }
            }
            assert(shrk_diag_rpt[loc_fi] + d_cnt == shrk_diag_rpt[loc_fi + 1]);
            for (idx_t p = S_mat.offd.row_ptr[loc_fi]; p < S_mat.offd.row_ptr[loc_fi + 1]; p ++) {
                if (S_mat.offd.col_idx[p] != -1) {
                    std::pair<int, idx_t> tmp;
                    tmp.first  = -1;
                    tmp.second = -1;
                    glb2pidloc.emplace(S_mat.col_map_offd[S_mat.offd.col_idx[p]], tmp);

                    shrk_offd_cid[shrk_offd_rpt[loc_fi] + o_cnt] = S_mat.col_map_offd[S_mat.offd.col_idx[p]];// 先改为全局序号暂存
                    shrk_offd_val[shrk_offd_rpt[loc_fi] + o_cnt] = S_mat.offd.vals   [p];
                    o_cnt ++;
                }
            }
            assert(shrk_offd_rpt[loc_fi] + o_cnt == shrk_offd_rpt[loc_fi + 1]);
        }
        delete [] S_mat.diag.row_ptr; S_mat.diag.row_ptr = shrk_diag_rpt;
        delete [] S_mat.diag.col_idx; S_mat.diag.col_idx = shrk_diag_cid;
        delete [] S_mat.diag.vals   ; S_mat.diag.vals    = shrk_diag_val;
        delete [] S_mat.offd.row_ptr; S_mat.offd.row_ptr = shrk_offd_rpt;
        delete [] S_mat.offd.col_idx; S_mat.offd.col_idx = shrk_offd_cid;
        delete [] S_mat.offd.vals   ; S_mat.offd.vals    = shrk_offd_val;
        S_mat.offd.ncols = glb2pidloc.size();
        // 置空
        S_mat.commpkg = nullptr;
        S_mat.col_map_offd = nullptr;
        // 会建立通信包和非对角的列映射
        S_mat.prepare_commpkg(glb2pidloc);
        // 再重新更新非对角部分的局部列号
        for (idx_t p = 0; p < S_mat.offd.nnz; p ++) {
            const idx_t glb_j = S_mat.offd.col_idx[p];
            assert(glb2pidloc.find(glb_j) != glb2pidloc.end());
            S_mat.offd.col_idx[p] = glb2pidloc[ glb_j ].second;
        }
    }
    delete [] my_dvs; my_dvs = nullptr;

#elif defined(TOPK)
    const float top_ratio = theta;
    idx_t * shrk_diag_rpt = new idx_t [my_num_F + 1]; shrk_diag_rpt[0] = 0;
    idx_t * shrk_offd_rpt = new idx_t [my_num_F + 1]; shrk_offd_rpt[0] = 0;
    for (idx_t i = 0; i < my_num_F; i++) {// 逐行计算Frobenius
        idx_t max_cnt = S_mat.diag.row_ptr[i + 1] - S_mat.diag.row_ptr[i]
                    +   S_mat.offd.row_ptr[i + 1] - S_mat.offd.row_ptr[i];
        idx_t glb_ids [max_cnt];
        setup_t vals  [max_cnt];
        idx_t cnt = 0;
        for (idx_t p = S_mat.diag.row_ptr[i]; p < S_mat.diag.row_ptr[i + 1]; p ++) {
            S_mat.diag.col_idx[p] = par_fine_mat->diag.col_idx[p];
            S_mat.diag.vals   [p] = calc_blk_norm(par_fine_mat->diag.vals + p * e_size);
            glb_ids[cnt] = S_mat.diag.col_idx[p] + _fbeg;
            vals   [cnt] = S_mat.diag.vals   [p];
            cnt ++;
        }
        for (idx_t p = S_mat.offd.row_ptr[i]; p < S_mat.offd.row_ptr[i + 1]; p ++) {
            S_mat.offd.col_idx[p] = par_fine_mat->offd.col_idx[p];
            S_mat.offd.vals   [p] = calc_blk_norm(par_fine_mat->offd.vals + p * e_size);
            glb_ids[cnt] = S_mat.col_map_offd[ S_mat.offd.col_idx[p] ];
            vals   [cnt] = S_mat.offd.vals   [p];
            cnt ++;
        }
        assert(cnt == max_cnt);
        sort_vals_with_indices<idx_t, setup_t>(max_cnt, vals, glb_ids);
        idx_t cut_off = max_cnt;
        // 先从后往前排除值为0的
        while (vals[cut_off - 1] == 0.0) cut_off --;
        // 如果剩下的仍比top_ratio多，则
        cut_off = MIN(cut_off, max_cnt * top_ratio);
        std::unordered_set<idx_t> remain;
        for (idx_t j = 0; j < cut_off; j++) remain.emplace(glb_ids[j]);
        idx_t dz_cnt = 0, oz_cnt = 0;
        for (idx_t p = S_mat.diag.row_ptr[i]; p < S_mat.diag.row_ptr[i + 1]; p ++) {
            const idx_t gid = S_mat.diag.col_idx[p] + _fbeg;
            if (remain.find(gid) == remain.end()) { S_mat.diag.col_idx[p] = -1; dz_cnt ++; }
        }
        for (idx_t p = S_mat.offd.row_ptr[i]; p < S_mat.offd.row_ptr[i + 1]; p ++) {
            const idx_t gid = S_mat.col_map_offd[ S_mat.offd.col_idx[p] ];
            if (remain.find(gid) == remain.end()) { S_mat.offd.col_idx[p] = -1; oz_cnt ++; }
        }
        shrk_diag_rpt[i + 1] = S_mat.diag.row_ptr[i + 1] - S_mat.diag.row_ptr[i] - dz_cnt;
        shrk_offd_rpt[i + 1] = S_mat.offd.row_ptr[i + 1] - S_mat.offd.row_ptr[i] - oz_cnt;
    }
    for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi ++) {
        shrk_diag_rpt[loc_fi + 1] += shrk_diag_rpt[loc_fi];
        shrk_offd_rpt[loc_fi + 1] += shrk_offd_rpt[loc_fi];
    }
    S_mat.diag.nnz = shrk_diag_rpt[my_num_F];
    S_mat.offd.nnz = shrk_offd_rpt[my_num_F];
    idx_t * shrk_diag_cid = new idx_t [S_mat.diag.nnz]; setup_t * shrk_diag_val = new setup_t [S_mat.diag.nnz];
    idx_t * shrk_offd_cid = new idx_t [S_mat.offd.nnz]; setup_t * shrk_offd_val = new setup_t [S_mat.offd.nnz];
    std::unordered_map<idx_t, std::pair<int, idx_t> > glb2pidloc;
    for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi ++) {
        idx_t d_cnt = 0, o_cnt = 0;
        for (idx_t p = S_mat.diag.row_ptr[loc_fi]; p < S_mat.diag.row_ptr[loc_fi + 1]; p ++) {
            if (S_mat.diag.col_idx[p] != -1) {
                shrk_diag_cid[shrk_diag_rpt[loc_fi] + d_cnt] = S_mat.diag.col_idx[p];
                shrk_diag_val[shrk_diag_rpt[loc_fi] + d_cnt] = S_mat.diag.vals   [p];
                d_cnt ++;
            }
        }
        assert(shrk_diag_rpt[loc_fi] + d_cnt == shrk_diag_rpt[loc_fi + 1]);
        for (idx_t p = S_mat.offd.row_ptr[loc_fi]; p < S_mat.offd.row_ptr[loc_fi + 1]; p ++) {
            if (S_mat.offd.col_idx[p] != -1) {
                std::pair<int, idx_t> tmp;
                tmp.first  = -1;
                tmp.second = -1;
                glb2pidloc.emplace(S_mat.col_map_offd[S_mat.offd.col_idx[p]], tmp);

                shrk_offd_cid[shrk_offd_rpt[loc_fi] + o_cnt] = S_mat.col_map_offd[S_mat.offd.col_idx[p]];// 先改为全局序号暂存
                shrk_offd_val[shrk_offd_rpt[loc_fi] + o_cnt] = S_mat.offd.vals   [p];
                o_cnt ++;
            }
        }
        assert(shrk_offd_rpt[loc_fi] + o_cnt == shrk_offd_rpt[loc_fi + 1]);
    }
    delete [] S_mat.diag.row_ptr; S_mat.diag.row_ptr = shrk_diag_rpt;
    delete [] S_mat.diag.col_idx; S_mat.diag.col_idx = shrk_diag_cid;
    delete [] S_mat.diag.vals   ; S_mat.diag.vals    = shrk_diag_val;
    delete [] S_mat.offd.row_ptr; S_mat.offd.row_ptr = shrk_offd_rpt;
    delete [] S_mat.offd.col_idx; S_mat.offd.col_idx = shrk_offd_cid;
    delete [] S_mat.offd.vals   ; S_mat.offd.vals    = shrk_offd_val;
    S_mat.offd.ncols = glb2pidloc.size();
    // 置空
    S_mat.commpkg = nullptr;
    S_mat.col_map_offd = nullptr;
    // 会建立通信包和非对角的列映射
    S_mat.prepare_commpkg(glb2pidloc);
    // 再重新更新非对角部分的局部列号
    for (idx_t p = 0; p < S_mat.offd.nnz; p ++) {
        const idx_t glb_j = S_mat.offd.col_idx[p];
        assert(glb2pidloc.find(glb_j) != glb2pidloc.end());
        S_mat.offd.col_idx[p] = glb2pidloc[ glb_j ].second;
    }

#elif defined(NORMAL)
    for (idx_t p = 0; p <  S_mat.diag.row_ptr[my_num_F]; p ++) {
        S_mat.diag.col_idx[p] = par_fine_mat->diag.col_idx[p];
        S_mat.diag.vals   [p] = calc_blk_norm(par_fine_mat->diag.vals + p * e_size);
    }
    for (idx_t p = 0; p <  S_mat.offd.row_ptr[my_num_F]; p ++) {
        S_mat.offd.col_idx[p] = par_fine_mat->offd.col_idx[p];
        S_mat.offd.vals   [p] = calc_blk_norm(par_fine_mat->offd.vals + p * e_size);
    }
    
#endif
    sort_and_find_bdr();

    assert(sorted);
    par_CSRMatrix<idx_t, setup_t, setup_t, 1> S_trans(comm);
    par_CSRMatrixTranspose(S_mat, S_trans, false);

    // // 打印强度矩阵出来看看
    // for (int p = 0; p < num_procs; p++) {
    //     if (my_pid == p) {
    //         printf("Proc %d : [%d, %d)\n", my_pid, _fbeg, _fend);
    //         for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi ++) {
    //             setup_t max_v = S_mat.diag.vals[S_mat.diag.row_ptr[loc_fi]];
    //             idx_t max_id = S_mat.diag.col_idx[S_mat.diag.row_ptr[loc_fi]] + _fbeg;
    //             if (S_mat.offd.row_ptr[loc_fi + 1] > S_mat.offd.row_ptr[loc_fi] &&
    //                 max_v < S_mat.offd.vals[S_mat.offd.row_ptr[loc_fi]] ) {
    //                 max_v = S_mat.offd.vals[S_mat.offd.row_ptr[loc_fi]];
    //                 max_id= S_mat.col_map_offd[ S_mat.offd.col_idx[S_mat.offd.row_ptr[loc_fi]] ];
    //             }
    //             if (max_id != loc_fi + _fbeg) printf("Warning !! ");
    //             for (idx_t p = S_mat.diag.row_ptr[loc_fi]; p < S_mat.diag.row_ptr[loc_fi + 1]; p ++) {
    //                 idx_t id = S_mat.diag.col_idx[p] + _fbeg;
    //                 setup_t v= S_mat.diag.vals   [p] / max_v;
    //                 printf("(%d, %.3f) ", id, v);
    //             }
    //             printf(" | ");
    //             for (idx_t p = S_mat.offd.row_ptr[loc_fi]; p < S_mat.offd.row_ptr[loc_fi + 1]; p ++) {
    //                 idx_t id = S_mat.col_map_offd[ S_mat.offd.col_idx[p] ];
    //                 setup_t v= S_mat.offd.vals   [p] / max_v;
    //                 printf("(%d, %.3f) ", id, v);
    //             }
    //             printf("\n");
    //         }
    //     }
    //     MPI_Barrier(comm);
    // } 

    const idx_t my_num_bdr_F = bdr_F.size();
    
    std::set<int> lower_pids;// 邻居里进程号小于自己的进程，注意这里是对称图
    std::set<int> upper_pids;// 邻居里进程号大于自己的进程，注意这里是对称图
    // int proc_dom[num_procs];
    // const int dom_size = num_procs / 32;
    // for (int p = 0; p < num_procs; p++) proc_dom[p] = p / dom_size;
    for (int r = 0; r < S_trans.commpkg->num_recvs; r++) {
        const int pid = S_trans.commpkg->recv_pids[r];
        // if (proc_dom[pid] != proc_dom[my_pid]) continue;
        if (pid < my_pid) lower_pids.emplace(pid);
        else              upper_pids.emplace(pid);
    }
    for (int s = 0; s < S_trans.commpkg->num_sends; s++) {
        const int pid = S_trans.commpkg->send_pids[s];
        // if (proc_dom[pid] != proc_dom[my_pid]) continue;
        if (pid < my_pid) lower_pids.emplace(pid);
        else              upper_pids.emplace(pid);
    }

    {// resolve 2-hop 
        const par_CSRCommPkg<idx_t> * commpkg_F = S_mat.commpkg;// par_fine_mat->commpkg;
        std::map<idx_t, std::unordered_set<int> > elmt2proc;
        for (int s = 0; s < commpkg_F->num_sends; s++) {
            const int pid = commpkg_F->send_pids[s];
            for (idx_t se = commpkg_F->send_map_starts[s]; se < commpkg_F->send_map_starts[s + 1]; se ++) {
                idx_t se_id = commpkg_F->send_map_elmts[se];
                if (elmt2proc.find(se_id) == elmt2proc.end()) elmt2proc.emplace(se_id, std::unordered_set<int>());
                elmt2proc[se_id].emplace(pid);
            }
        }
        // 从中去除只有1个进程的集合
        std::map<int, idx_t> proc_cnts;// 记录给每个进程需要发送的缓冲区长度
        for (auto it = elmt2proc.begin(); it != elmt2proc.end(); ) {
            if (it->second.size() == 1) { elmt2proc.erase(it ++); }// 满足删除条件，删除当前结点，并指向下面一个结点
            else {
                for (auto jt = it->second.begin(); jt != it->second.end(); jt ++) {
                    if (proc_cnts.find(*jt) == proc_cnts.end()) proc_cnts.emplace(*jt, 0);
                    proc_cnts[*jt] += (it->second.size() - 1);
                }
                it ++;// 条件不满足，指向下面一个结点
            }
        }
        idx_t buf_len = 0;
        for (auto it = proc_cnts.begin(); it != proc_cnts.end(); it ++) buf_len += it->second;

        idx_t * send_buf_0 = new idx_t [commpkg_F->num_sends];
        int   * send_buf_1 = new int   [buf_len];
        idx_t buf1_ptr = 0;
        std::vector<MPI_Request> send_reqs;
        for (idx_t s = 0; s < commpkg_F->num_sends; s++) {
            const int pid = commpkg_F->send_pids[s];
            if (proc_cnts.find(pid) == proc_cnts.end()) send_buf_0[s] = 0;
            else                                        send_buf_0[s] = proc_cnts[pid];
            MPI_Request req_num, req_dat;
            MPI_Isend(& send_buf_0[s], 1, mpi_idx_type, commpkg_F->send_pids[s], 89, commpkg_F->comm, & req_num);
            send_reqs.emplace_back(req_num);
            if (send_buf_0[s] > 0) {
                idx_t my_ptr = buf1_ptr;
                for (auto it = elmt2proc.begin(); it != elmt2proc.end(); it ++) {
                    if (it->second.find(pid) != it->second.end()) {
                        for (auto jt = it->second.begin(); jt != it->second.end(); jt ++) {
                            if (*jt != pid) send_buf_1[my_ptr ++] = *jt;
                        }
                    }
                }
                assert(my_ptr == buf1_ptr + send_buf_0[s]);
                MPI_Isend(send_buf_1 + buf1_ptr, send_buf_0[s], MPI_INT, commpkg_F->send_pids[s], 90, commpkg_F->comm, & req_dat);
                send_reqs.emplace_back(req_dat);
            }
            buf1_ptr += send_buf_0[s];
        }
        assert(buf1_ptr == buf_len);

        for (idx_t r = 0; r < commpkg_F->num_recvs; r++) {
            idx_t cnt = 0;
            MPI_Recv(& cnt, 1, mpi_idx_type, commpkg_F->recv_pids[r], 89, commpkg_F->comm, MPI_STATUS_IGNORE);
            assert(cnt >= 0);
            if (cnt == 0) continue;
            int * buf = new int [cnt];
            MPI_Recv(buf, cnt, MPI_INT, commpkg_F->recv_pids[r], 90, commpkg_F->comm, MPI_STATUS_IGNORE);
            for (idx_t i = 0; i < cnt; i++) {
                int pid = buf[i]; assert(pid != my_pid);
                if (pid < my_pid) lower_pids.emplace(pid);
                else              upper_pids.emplace(pid);
            }
            delete [] buf;
        }
        MPI_Waitall(send_reqs.size(), send_reqs.data(), MPI_STATUSES_IGNORE);
        delete [] send_buf_0; delete [] send_buf_1;
#ifdef DEBUG
    int print_procs[] = {131, 240, 241};
    for (int p = 0; p < sizeof(print_procs)/sizeof(int); p ++) {
        if (print_procs[p] == my_pid) {
            printf("Proc %d lower : ", my_pid);
            for (auto it = lower_pids.begin(); it != lower_pids.end(); it ++) printf("%d ", *it);
            printf("\nProc %d upper : ", my_pid);
            for (auto it = upper_pids.begin(); it != upper_pids.end(); it ++) printf("%d ", *it);
            printf("\nProc %d recv_pids : ", my_pid);
            for (int r = 0; r < commpkg_F->num_recvs; r++) printf("%d ", commpkg_F->recv_pids[r]);
            printf("\nProc %d send_pids : ", my_pid);
            for (int s = 0; s < commpkg_F->num_sends; s++) printf("%d ", commpkg_F->send_pids[s]);
            printf("\n");
        }
        MPI_Barrier(MPI_COMM_WORLD);
    }
#endif
    }
    
    // {// 统计一下邻居进程数、边界点数
    //     long long _loc[3];
    //     _loc[0] = bdr_F.size();
    //     _loc[1] = lower_pids.size();
    //     _loc[2] = upper_pids.size();
    //     long long   _glb_min[sizeof(_loc)/sizeof(long long)],
    //                 _glb_max[sizeof(_loc)/sizeof(long long)],
    //                 _glb_avg[sizeof(_loc)/sizeof(long long)];
    //     MPI_Reduce(_loc, _glb_min, 3, MPI_LONG_LONG, MPI_MIN, 0, comm);
    //     MPI_Reduce(_loc, _glb_avg, 3, MPI_LONG_LONG, MPI_SUM, 0, comm);
    //     MPI_Reduce(_loc, _glb_max, 3, MPI_LONG_LONG, MPI_MAX, 0, comm);
    //     for (size_t i = 0; i < sizeof(_loc)/sizeof(long long); i++) _glb_avg[i] /= num_procs;
    //     if (my_pid == 0) {
    //         printf("Min %lld %lld %lld\n", _glb_min[0], _glb_min[1], _glb_min[2]);
    //         printf("Avg %lld %lld %lld\n", _glb_avg[0], _glb_avg[1], _glb_avg[2]);
    //         printf("Max %lld %lld %lld\n", _glb_max[0], _glb_max[1], _glb_max[2]);
    //     }
    // }
    
    // -2: aggregated by other procs, -1: unaggregated, >= 0: aggregated by myself with this coar local ID
    idx_t * loc_map_f2c = new idx_t [my_num_F];
    coar_wells.clear();// 为这一轮粗化清空
    constexpr idx_t UNAGGREGATED = -1, AGGREGATED_BY_OTHERS = -2;
    // 初始化
    for (idx_t fi = 0; fi < my_num_F; fi++)
        loc_map_f2c[fi] = UNAGGREGATED;

    std::unordered_map<idx_t, idx_t> col_map_glb_FT;
    for (idx_t lojT = 0; lojT < S_trans.offd.ncols; lojT ++) {// lojT: local offd j of S_trans
        const idx_t glb_fjT = S_trans.col_map_offd[lojT];
        col_map_glb_FT.emplace(glb_fjT, lojT);
    }

    // double t2 = wall_time();

    idx_t phase1_C = 0;
    {// 先聚合自己处于进程边界的点
        struct STATUS {
            idx_t id = -1;// 全局序号
            idx_t agg_id = -1;// 与其聚合的全局序号
        };
        std::vector<MPI_Request> send_reqs;
        
        idx_t send_cnt = 0;
        // 发送缓冲区，记录处于本进程的边界点的状态：-1: unaggregated, 其它: 与其聚合的邻居的全局序号
        std::vector<STATUS> send_buf;
        
        // 记录处于别的进程的边界点的状态，-1: unaggregated, 其它：与其聚合的邻居的全局序号
        std::unordered_map<idx_t, idx_t> ngb_F_status;
        for (idx_t j = 0; j < S_mat.offd.ncols; j ++) {// 初始化为unaggregated
            ngb_F_status.emplace(S_mat.col_map_offd[j], UNAGGREGATED);
            // printf("Record %d => %d\n", S_mat.col_map_offd[j], ngb_F_status[S_mat.col_map_offd[j]]);
        }

        for (int work_pid = 0; work_pid < num_procs; work_pid ++) {// 边界点的聚合按进程编号从小到大进行
            if (work_pid == my_pid) {
                for (idx_t i = 0; i < my_num_bdr_F; i ++) {
                    const idx_t loc_fi = bdr_F[i];
                    STATUS stat;
                    stat.id = _fbeg + loc_fi;

                    if (loc_map_f2c[loc_fi] <= AGGREGATED_BY_OTHERS) {// 之前已经被进程号比自己小的进程聚合过了
                        stat.agg_id = AGGREGATED_BY_OTHERS;
                    }
                    else if (loc_map_f2c[loc_fi] > 0) {// 被之前的某个自己进程负责的边界点聚合过了
                        // 不用管，该agg_id一定在之前的某个stat.id中
                    }
                    else {
                        // 找出diag内最强耦合的邻居
                        idx_t diag_candi = -1; setup_t diag_v = 0.0;// idx_t dc_ptr = -1;
                        for (idx_t p = S_mat.diag.row_ptr[loc_fi]; p < S_mat.diag.row_ptr[loc_fi + 1]; p ++) {// 从大到小考察邻居
                            const idx_t loc_fj = S_mat.diag.col_idx[p];
                            const idx_t glb_fj = loc_fj + _fbeg;
                            // dc_ptr = p - S_mat.diag.row_ptr[loc_fi];
                            if (loc_fj != loc_fi && // 该邻居不能是自己
                                iso_wells.find(glb_fj) == iso_wells.end() && // 也不能是井
                                loc_map_f2c[loc_fj] == UNAGGREGATED) // 也尚未被聚合
                            {
                                diag_candi = glb_fj;
                                diag_v     = S_mat.diag.vals[p];
                                break;
                            }
                        }
                        // 找出offd内最强耦合的邻居
                        idx_t offd_candi = -1; setup_t offd_v = 0.0;// idx_t oc_ptr = -1;
                        for (idx_t p = S_mat.offd.row_ptr[loc_fi]; p < S_mat.offd.row_ptr[loc_fi + 1]; p ++) {// 从大到小考察邻居
                            const idx_t loc_fj = S_mat.offd.col_idx[p];
                            const idx_t glb_fj = S_mat.col_map_offd[loc_fj];
                            // oc_ptr = p - S_mat.offd.row_ptr[loc_fi];
                            if (ngb_F_status[glb_fj] == UNAGGREGATED && // 该邻居尚未被（自己或更小的进程）聚合
                                iso_wells.find(glb_fj) == iso_wells.end() )// 也不能是井
                            {
                                offd_candi = glb_fj;
                                offd_v     = S_mat.offd.vals[p];
                                break;
                            }
                        }
                        idx_t final_candi = -1;// 最终确定的聚合邻居
                        if (diag_candi != -1 && offd_candi != -1) {
                            if (diag_v > offd_v) final_candi = diag_candi;// 选值大的
                            else                 final_candi = offd_candi;
                        }
                        else { final_candi = MAX(diag_candi, offd_candi); }

                        if (final_candi == -1) {// 找不到合适的邻居进行聚合
                            // printf("Warning !! Proc %d loc %d glb %d at bdr got not aggre! num_diag %d (%d) num_offd %d (%d)\n",
                            //     my_pid, loc_fi, stat.id,    S_mat.diag.row_ptr[loc_fi + 1] - S_mat.diag.row_ptr[loc_fi], diag_candi,
                            //                                 S_mat.offd.row_ptr[loc_fi + 1] - S_mat.offd.row_ptr[loc_fi], offd_candi);
                            // for (idx_t p = S_mat.diag.row_ptr[loc_fi]; p < S_mat.diag.row_ptr[loc_fi + 1]; p ++) {
                            //     const idx_t loc_fj = S_mat.diag.col_idx[p];
                            //     if (loc_fj != loc_fi)
                            //         printf("glb %d => %d\n", S_mat.diag.col_idx[p] + _fbeg, loc_map_f2c[loc_fj]);
                            // }
                            // for (idx_t p = S_mat.offd.row_ptr[loc_fi]; p < S_mat.offd.row_ptr[loc_fi + 1]; p ++) {
                            //     const idx_t loc_fj = S_mat.offd.col_idx[p];
                            //     const idx_t glb_fj = S_mat.col_map_offd[loc_fj];
                            //     printf("glb %d => %d\n", S_mat.col_map_offd[loc_fj], ngb_F_status[glb_fj]);
                            // }
                        }
                        else {
                            // printf("Proc %d Glb %d pick %d\n", my_pid, stat.id, final_candi);
                            if (final_candi == diag_candi) { assert(diag_candi != -1); // 聚合的邻居是自己本进程内的
                                const idx_t loc_fj = diag_candi - _fbeg;
                                loc_map_f2c[loc_fi] = my_num_F + phase1_C;
                                loc_map_f2c[loc_fj] = my_num_F + phase1_C;
                            } else { assert(offd_candi != -1); // 聚合的邻居位于别的进程
                                assert(ngb_F_status.find(offd_candi) != ngb_F_status.end() && ngb_F_status[offd_candi] == UNAGGREGATED);
                                loc_map_f2c[loc_fi]      = my_num_F + phase1_C;
                                ngb_F_status[offd_candi] = my_num_F + phase1_C;
                            }
                            phase1_C ++;

                            stat.agg_id = final_candi;
                        }   
                    }
                    send_buf.emplace_back(stat);
                }// boundary F loop
                assert(send_buf.size() == bdr_F.size());
#ifdef DEBUG
                printf("Proc %d Finish bdr\n", my_pid);
                for (int loc_fi = 0; loc_fi < my_num_F; loc_fi ++) {
                    printf("   glb %d => %d\n", loc_fi + _fbeg, loc_map_f2c[loc_fi]);
                }
#endif
                send_cnt = send_buf.size();
                // 向进程比自己大且为邻居的进程发送消息更新它们负责范围内的边界点（是否已被聚合）
                for (auto it = upper_pids.begin(); it != upper_pids.end(); it ++) {
                    const int dst_pid = *it;
                    MPI_Request req_cnt, req_dat;
                    MPI_Isend(& send_cnt, 1, mpi_idx_type, dst_pid, 91, comm, & req_cnt);
                    MPI_Isend(send_buf.data(), send_cnt * sizeof(STATUS), MPI_BYTE, dst_pid, 92, comm, & req_dat);
                    send_reqs.emplace_back(req_cnt);
                    send_reqs.emplace_back(req_dat);
                }
            }// my job
            else if (work_pid < my_pid && lower_pids.find(work_pid) != lower_pids.end()) {// 该小进程号的进程是我的邻居
                idx_t cnt = 0;
                MPI_Recv(& cnt, 1, mpi_idx_type, work_pid, 91, comm, MPI_STATUS_IGNORE);
                STATUS * recv_buf = new STATUS [cnt];// 收取work_pid这个进程对它自身的边界点的记录
                MPI_Recv(recv_buf, cnt * sizeof(STATUS), MPI_BYTE, work_pid, 92, comm, MPI_STATUS_IGNORE);
                
                // update status
                for (idx_t k = 0; k < cnt; k++) {
                    const idx_t glb_id = recv_buf[k].id;// 属于work_pid进程负责的点（全局序号）
                    const idx_t agg_id = recv_buf[k].agg_id;
                    // if (agg_id == 54091896) {
                    //     printf("Proc %d from Proc %d : glb_id %d agg_id %d\n", my_pid, work_pid, glb_id, agg_id);
                    // }
                    if (agg_id != -1) {// work_pid这个进程将它范围内的glb_id这个边界点聚合了，
                        // 如果glb_id这个边界点属于与我邻接，则需要将它的状态更新
                        if (ngb_F_status.find(glb_id) != ngb_F_status.end()) {
                            if (ngb_F_status[glb_id] != UNAGGREGATED) {// 记录显示该glb_id的点之前已经被聚合过
                                assert(agg_id == AGGREGATED_BY_OTHERS);
                            } else {// 记录仍未未聚合
                                ngb_F_status[glb_id] = agg_id;
                            }
                        }
                        // 进一步检查我是否有关于agg_id的信息记录
                        // work_pid这个进程将属于我负责的agg_id这个点与它负责的glb_id这个点聚合了，故需要更新我的局部映射信息
                        if (_fbeg <= agg_id && agg_id < _fend) {
                            // if (col_map_glb_FT.find(glb_id) == col_map_glb_FT.end()) {
                            //     printf("Proc %d [%d,%d) Error !! from Proc %d glb_id %d agg_id %d cannot find in FT\n", my_pid,
                            //         _fbeg, _fend, work_pid, glb_id, agg_id);
                            // }
                            assert(col_map_glb_FT.find(glb_id) != col_map_glb_FT.end());// glb_id这个点应属于与我相邻的边界点
                            loc_map_f2c[agg_id - _fbeg] = AGGREGATED_BY_OTHERS - glb_id;// 注意减掉
                        }
                        // work_pid这个进程将属于我的另一邻居进程（可能也是work_pid）的agg_id这个点与它负责的glb_id这个点聚合了，
                        // 故需要更新我对别的进程的边界点的记录
                        else if (ngb_F_status.find(agg_id) != ngb_F_status.end()) {
                            // if (ngb_F_status[agg_id] != UNAGGREGATED) {
                            //     printf("Proc %d Error !! from Proc %d glb_id %d agg_id %d ngb_F_status %d\n", my_pid, work_pid,
                            //         glb_id, agg_id, ngb_F_status[agg_id]);
                            // }
                            assert(ngb_F_status[agg_id] == UNAGGREGATED);
                            ngb_F_status[agg_id] = glb_id;
                        }
                    }
                }
                delete [] recv_buf;
            }
            else if (work_pid > my_pid && upper_pids.find(work_pid) != upper_pids.end()) {// 该大进程号的进程是我的邻居
                // nothing to do
            }
        }// pid loop

        MPI_Waitall(send_reqs.size(), send_reqs.data(), MPI_STATUSES_IGNORE);
    }

    // double t3 = wall_time();

    idx_t my_num_C = 0;
    std::unordered_map<idx_t, idx_t> map_phase1C_coarID;
    {// 再聚合内部点
        for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi++) {
            if (loc_map_f2c[loc_fi] != UNAGGREGATED) {// 之前已经聚合过了
                if (loc_map_f2c[loc_fi] >= my_num_F) {// 是第一阶段被自己聚合的边界点
                    if (map_phase1C_coarID.find(loc_map_f2c[loc_fi]) == map_phase1C_coarID.end()) {
                        map_phase1C_coarID.emplace(loc_map_f2c[loc_fi], my_num_C);
                        loc_map_f2c[loc_fi] = my_num_C;
                        my_num_C ++;
                    } else {
                        loc_map_f2c[loc_fi] = map_phase1C_coarID[loc_map_f2c[loc_fi]];
                    }
                }// 否则 [0, my_num_F)的是刚才被自己聚合的，<= UNAGGREGATED(-1)即只能是被其它进程聚合的
                continue;
            }
            if (iso_wells.find(loc_fi + _fbeg) != iso_wells.end()){// 井不与别人聚合，单独成一个粗点
                assert(loc_map_f2c[loc_fi] == UNAGGREGATED);// 之前不能被聚合过
                coar_wells.emplace(my_num_C);
                loc_map_f2c[loc_fi] = my_num_C;
                my_num_C ++;
                continue;
            }

            // 从大到小考察邻居
            for (idx_t p = S_mat.diag.row_ptr[loc_fi]; p < S_mat.diag.row_ptr[loc_fi + 1]; p++) {
                const idx_t loc_fj = S_mat.diag.col_idx[p];
                if (loc_fj != loc_fi && // 该邻居不能是自己
                    iso_wells.find(loc_fj + _fbeg) == iso_wells.end() && // 也不能是井
                    loc_map_f2c[loc_fj] == UNAGGREGATED) // 也尚未被聚合
                {
                    loc_map_f2c[loc_fi] = my_num_C;
                    loc_map_f2c[loc_fj] = my_num_C;
                    // printf("%d %d => %d\n", loc_fi, loc_fj, my_num_C);
                    my_num_C ++;
                    break;
                }
            }
        }// fi loop
    }
#ifdef DEBUG
    for (int p = 0; p < num_procs; p++) {
        if (my_pid == p) {
            printf("Proc %d my_num_bdr_F %d\n", my_pid, my_num_bdr_F);
            for (int j = 0; j < my_num_bdr_F; j++) {
                const idx_t loc_fi = bdr_F[j];
                printf("   glb %d => %d\n", loc_fi + _fbeg, loc_map_f2c[loc_fi]);
            }
        }
        MPI_Barrier(comm);
    }
#endif
#if 0
    for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi++) {
        // 自成一点，会使粗层网格数偏多，且效果不一定更好
        if (loc_map_f2c[loc_fi] == UNAGGREGATED) {
            assert(iso_wells.find(loc_fi + _fbeg) == iso_wells.end());
            loc_map_f2c[loc_fi] = my_num_C ++;
        }
    }
#endif
    // double t4 = wall_time();

    idx_t * C_rows_partits = new idx_t [num_procs + 1];
    C_rows_partits[0] = 0;
    MPI_Allgather(& my_num_C, 1, mpi_idx_type, C_rows_partits + 1, 1, mpi_idx_type, comm);
    for (int p = 0; p < num_procs; p++)
        C_rows_partits[p + 1] += C_rows_partits[p];

    const idx_t _cbeg = C_rows_partits[my_pid], _cend = C_rows_partits[my_pid + 1];
    idx_t * glb_map_f2c = new idx_t [my_num_F];// loc_f => glb_c
    for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi ++) {
        if (loc_map_f2c[loc_fi] >= 0)
            glb_map_f2c[loc_fi] = loc_map_f2c[loc_fi] + _cbeg;// 累加全局偏移
        else
            glb_map_f2c[loc_fi] = -1;// 暂时置为-1
    }
    {// 告知邻居进程，细网格上我的边界点的全局序号
        idx_t * send_buf = new idx_t [S_trans.commpkg->send_map_starts[S_trans.commpkg->num_sends]];
        idx_t * recv_buf = new idx_t [S_trans.offd.ncols];
        for (idx_t s = 0; s < S_trans.commpkg->num_sends; s++) {
            idx_t       * buf    = send_buf                        + S_trans.commpkg->send_map_starts[s];
            const idx_t * my_ids = S_trans.commpkg->send_map_elmts + S_trans.commpkg->send_map_starts[s];
            const idx_t num = S_trans.commpkg->send_map_starts[s + 1] - S_trans.commpkg->send_map_starts[s];
            for (idx_t k = 0; k < num; k++) {
                buf[k] = glb_map_f2c[my_ids[k]];
            }
            MPI_Isend(buf, num, mpi_idx_type, S_trans.commpkg->send_pids[s], 113, comm, S_trans.commpkg->send_reqs + s);
        }
        for (idx_t r = 0; r < S_trans.commpkg->num_recvs; r++) {
            idx_t * buf = recv_buf + S_trans.commpkg->recv_vec_starts[r];
            const idx_t num = S_trans.commpkg->recv_vec_starts[r + 1] - S_trans.commpkg->recv_vec_starts[r];
            MPI_Irecv(buf, num, mpi_idx_type, S_trans.commpkg->recv_pids[r], 113, comm, S_trans.commpkg->recv_reqs + r);
        }
        MPI_Waitall(S_trans.commpkg->num_recvs, S_trans.commpkg->recv_reqs, MPI_STATUSES_IGNORE);

        for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi++) {
            if (loc_map_f2c[loc_fi] <= AGGREGATED_BY_OTHERS) { assert(glb_map_f2c[loc_fi] == -1);
                const idx_t glb_fj = - (loc_map_f2c[loc_fi] - AGGREGATED_BY_OTHERS);
                assert(col_map_glb_FT.find(glb_fj) != col_map_glb_FT.end());
                const idx_t loc_fj = col_map_glb_FT[glb_fj];
                const idx_t glb_cj = recv_buf[ loc_fj ];
                glb_map_f2c[loc_fi] = glb_cj;// 更新全局映射
            }
            // 此时先不管 loc_map_f2c[loc_fi] == UNAGGREGATED 即尚未聚合的点
        }
        MPI_Waitall(S_trans.commpkg->num_sends, S_trans.commpkg->send_reqs, MPI_STATUSES_IGNORE);
        delete [] recv_buf; delete [] send_buf;
    }
    {// 处理尚未聚合的点
        idx_t * send_buf = new idx_t [S_mat.commpkg->send_map_starts[S_mat.commpkg->num_sends]];
        idx_t * recv_buf = new idx_t [S_mat.offd.ncols];
        for (idx_t s = 0; s < S_mat.commpkg->num_sends; s++) {
            idx_t       * buf    = send_buf                        + S_mat.commpkg->send_map_starts[s];
            const idx_t * my_ids = S_mat.commpkg->send_map_elmts + S_mat.commpkg->send_map_starts[s];
            const idx_t num = S_mat.commpkg->send_map_starts[s + 1] - S_mat.commpkg->send_map_starts[s];
            for (idx_t k = 0; k < num; k++) {
                buf[k] = glb_map_f2c[my_ids[k]];
            }
            MPI_Isend(buf, num, mpi_idx_type, S_mat.commpkg->send_pids[s], 114, comm, S_mat.commpkg->send_reqs + s);
        }
        for (idx_t r = 0; r < S_mat.commpkg->num_recvs; r++) {
            idx_t * buf = recv_buf + S_mat.commpkg->recv_vec_starts[r];
            const idx_t num = S_mat.commpkg->recv_vec_starts[r + 1] - S_mat.commpkg->recv_vec_starts[r];
            MPI_Irecv(buf, num, mpi_idx_type, S_mat.commpkg->recv_pids[r], 114, comm, S_mat.commpkg->recv_reqs + r);
        }
        MPI_Waitall(S_mat.commpkg->num_recvs, S_mat.commpkg->recv_reqs, MPI_STATUSES_IGNORE);

        for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi++) {
            if (loc_map_f2c[loc_fi] == UNAGGREGATED) {// 该点仍未被聚合
                // 从所有邻居中选最大的进行聚合：注意包含落在别的进程上的邻居
                // 找出diag内最强耦合的邻居
                idx_t diag_candi = -1; setup_t diag_v = 0.0;// idx_t dc_ptr = -1;
                for (idx_t p = S_mat.diag.row_ptr[loc_fi]; p < S_mat.diag.row_ptr[loc_fi + 1]; p ++) {// 从大到小考察邻居
                    const idx_t loc_fj = S_mat.diag.col_idx[p];
                    const idx_t glb_fj = loc_fj + _fbeg;
                    if (loc_fj != loc_fi && // 该邻居不能是自己
                        iso_wells.find(glb_fj) == iso_wells.end() ) // 也不能是井
                    {
                        assert(glb_map_f2c[loc_fj] != -1);
                        diag_candi = glb_map_f2c[loc_fj];
                        diag_v     = S_mat.diag.vals[p];
                        break;
                    }
                }
                // 找出offd内最强耦合的邻居
                idx_t offd_candi = -1; setup_t offd_v = 0.0;// idx_t oc_ptr = -1;
                for (idx_t p = S_mat.offd.row_ptr[loc_fi]; p < S_mat.offd.row_ptr[loc_fi + 1]; p ++) {// 从大到小考察邻居
                    const idx_t loc_fj = S_mat.offd.col_idx[p];
                    const idx_t glb_fj = S_mat.col_map_offd[loc_fj];
                    if (iso_wells.find(glb_fj) == iso_wells.end() )// 不能是井
                    {
                        assert(recv_buf[loc_fj] != -1);// 反证法可说明不应存在尚未被聚合的处于别的进程上的邻居
                        offd_candi = recv_buf[loc_fj];
                        offd_v     = S_mat.offd.vals[p];
                        break;
                    }
                }
                idx_t final_candi = -1;// 最终确定的聚合邻居
                if (diag_candi != -1 && offd_candi != -1) {
                    if (diag_v > offd_v) final_candi = diag_candi;// 选值大的
                    else                 final_candi = offd_candi;
                }
                else { final_candi = MAX(diag_candi, offd_candi); }
                // if (loc_fi + _fbeg == 405) printf("Proc %d glb_i %d diag_v %.2f offd_v %.2f\n", my_pid, loc_fi + _fbeg, diag_v, offd_v);
                glb_map_f2c[loc_fi] = final_candi;
                if (final_candi == -1) {// check: 此时还尚未被聚合的，一定是只有一个邻居（即自己的）
                    assert(S_mat.diag.row_ptr[loc_fi + 1] - S_mat.diag.row_ptr[loc_fi] == 1);
                    assert(S_mat.offd.row_ptr[loc_fi + 1] - S_mat.offd.row_ptr[loc_fi] == 0);
                    assert(S_mat.diag.col_idx[S_mat.diag.row_ptr[loc_fi]] == loc_fi);
                }
            }
        }
        MPI_Waitall(S_mat.commpkg->num_sends, S_mat.commpkg->send_reqs, MPI_STATUSES_IGNORE);
        delete [] recv_buf; delete [] send_buf;
    }// 此时在glb_map_f2c[]中仍为-1的点将不被粗化进入下一层粗网格

#ifdef DEBUG
    for (int p = 0; p < num_procs; p++) {
        if (my_pid == p) {
            printf("Proc %d num_F %d num_C %d\n", my_pid, my_num_F, my_num_C);
            for (int loc_fi = 0; loc_fi < my_num_F; loc_fi++) {
                const idx_t glb_fi = loc_fi + _fbeg;
                printf("  Glb %d locMap %d glbMap %d\n", glb_fi, loc_map_f2c[loc_fi], glb_map_f2c[loc_fi]);
            }
        }
        MPI_Barrier(comm);
    }
#endif
    // 构造限制矩阵
    R_mat = new par_CSRMatrix<idx_t, setup_t, setup_t, 1>(comm, C_rows_partits, par_fine_mat->rows_partition);
    P_mat = new par_CSRMatrix<idx_t, setup_t, setup_t, 1>(comm, par_fine_mat->cols_partition, C_rows_partits);

    {// 构造插值矩阵
        P_mat->diag.nrows = my_num_F; P_mat->offd.nrows = my_num_F;
        P_mat->diag.row_ptr = new idx_t [my_num_F + 1]; P_mat->diag.row_ptr[0] = 0;
        P_mat->offd.row_ptr = new idx_t [my_num_F + 1]; P_mat->offd.row_ptr[0] = 0;
        std::unordered_map<idx_t, std::pair<int, idx_t> > Poffd_g2l;// 记录插值矩阵的非对角部分的全局序到局部序的映射
        for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi ++) {
            const idx_t glb_ci = glb_map_f2c[loc_fi];
            if (_cbeg <= glb_ci && glb_ci < _cend) {// 粗点属于本进程负责计算
                P_mat->diag.row_ptr[loc_fi + 1] = 1;
                P_mat->offd.row_ptr[loc_fi + 1] = 0;
            } else if (glb_ci >= 0) {// 粗点不属于本进程计算
                P_mat->diag.row_ptr[loc_fi + 1] = 0;
                P_mat->offd.row_ptr[loc_fi + 1] = 1;
                
                std::pair<int, idx_t> tmp(-1, -1);
                Poffd_g2l.emplace(glb_ci, tmp);
            } else { assert(glb_ci == -1);// 粗点不进入下一层
                P_mat->diag.row_ptr[loc_fi + 1] = 0;
                P_mat->offd.row_ptr[loc_fi + 1] = 0;
            }
        }
        for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi ++) {// prefix sum
            P_mat->diag.row_ptr[loc_fi + 1] += P_mat->diag.row_ptr[loc_fi];
            P_mat->offd.row_ptr[loc_fi + 1] += P_mat->offd.row_ptr[loc_fi];
        }
        P_mat->prepare_commpkg(Poffd_g2l);

        // 确定列数和非零元数
        P_mat->diag.ncols = my_num_C;                      P_mat->offd.ncols = Poffd_g2l.size();
        P_mat->diag.nnz   = P_mat->diag.row_ptr[my_num_F]; P_mat->offd.nnz   = P_mat->offd.row_ptr[my_num_F];
        // 分配列序号和值的内存
        P_mat->diag.col_idx = new idx_t  [P_mat->diag.nnz]; P_mat->offd.col_idx = new idx_t  [P_mat->offd.nnz];
        P_mat->diag.vals    = new setup_t[P_mat->diag.nnz]; P_mat->offd.vals    = new setup_t[P_mat->offd.nnz];

        for (idx_t loc_fi = 0; loc_fi < my_num_F; loc_fi ++) {
            const idx_t glb_ci = glb_map_f2c[loc_fi];
            // const idx_t dcnt = 0, ocnt = 0;
            const idx_t doff = P_mat->diag.row_ptr[loc_fi],
                        ooff = P_mat->offd.row_ptr[loc_fi];
            if (_cbeg <= glb_ci && glb_ci < _cend) {// 粗点属于本进程负责计算
                P_mat->diag.col_idx[doff] = glb_ci - _cbeg;
                P_mat->diag.vals   [doff] = 1.0;
            } else if (glb_ci >= 0) {// 粗点不属于本进程计算
                assert(Poffd_g2l.find(glb_ci) != Poffd_g2l.end());
                P_mat->offd.col_idx[ooff] = Poffd_g2l[glb_ci].second;
                P_mat->offd.vals   [ooff] = 1.0;
            } else { assert(glb_ci == -1);// 粗点不进入下一层
                // nothing to do
            }
        }
#ifdef DEBUG
        for (int p = 0; p < num_procs; p++) {
            if (my_pid == p) {
                printf("Proc %d Prlg diag %d %d %d offd %d %d %d\n", my_pid, P_mat->diag.nrows, P_mat->diag.ncols, P_mat->diag.nnz,
                    P_mat->offd.nrows, P_mat->offd.ncols, P_mat->offd.nnz);
                for (idx_t i = 0; i < P_mat->diag.nrows; i++) {
                    printf("  %d : ", P_mat->beg_row + i);
                    for (idx_t p = P_mat->diag.row_ptr[i]; p < P_mat->diag.row_ptr[i + 1]; p++)
                        printf("%d ", P_mat->beg_col + P_mat->diag.col_idx[p]);
                    printf(" | ");
                    for (idx_t p = P_mat->offd.row_ptr[i]; p < P_mat->offd.row_ptr[i + 1]; p++)
                        printf("%d ", P_mat->col_map_offd[P_mat->offd.col_idx[p]]);
                    printf("\n");
                }
            }
            MPI_Barrier(comm);
        }
#endif
    }
    par_CSRMatrixTranspose(* P_mat, * R_mat);
    {// 需要修改限制矩阵的值
        for (idx_t i = 0; i < R_mat->diag.nrows; i ++) {
            const idx_t num = R_mat->diag.row_ptr[i + 1] - R_mat->diag.row_ptr[i]
                            + R_mat->offd.row_ptr[i + 1] - R_mat->offd.row_ptr[i];
            const setup_t v = 1.0 / (setup_t) num;
            for (idx_t p = R_mat->diag.row_ptr[i]; p < R_mat->diag.row_ptr[i + 1]; p ++) {
                assert(R_mat->diag.vals[p] == 1.0);
                R_mat->diag.vals[p] = v;
            }
            for (idx_t p = R_mat->offd.row_ptr[i]; p < R_mat->offd.row_ptr[i + 1]; p ++) {
                assert(R_mat->offd.vals[p] == 1.0);
                R_mat->offd.vals[p] = v;
            }
        }
#ifdef DEBUG
        for (int p = 0; p < num_procs; p++) {
            if (my_pid == p) {
                printf("Proc %d Rstr diag %d %d %d offd %d %d %d\n", my_pid, R_mat->diag.nrows, R_mat->diag.ncols, R_mat->diag.nnz,
                    R_mat->offd.nrows, R_mat->offd.ncols, R_mat->offd.nnz);
                for (idx_t i = 0; i < R_mat->diag.nrows; i++) {
                    printf("  %d : ", R_mat->beg_row + i);
                    for (idx_t p = R_mat->diag.row_ptr[i]; p < R_mat->diag.row_ptr[i + 1]; p++)
                        printf("%d ", R_mat->beg_col + R_mat->diag.col_idx[p]);
                    printf(" | ");
                    for (idx_t p = R_mat->offd.row_ptr[i]; p < R_mat->offd.row_ptr[i + 1]; p++)
                        printf("%d ", R_mat->col_map_offd[R_mat->offd.col_idx[p]]);
                    printf("\n");
                }
            }
            MPI_Barrier(comm);
        }
#endif
    }

// #ifdef DEBUG
    const idx_t my_num_I = my_num_F - P_mat->diag.nnz - P_mat->offd.nnz;// 独立点个数
    idx_t tot_num_I = 0;
    MPI_Reduce(& my_num_I, & tot_num_I, 1, mpi_idx_type, MPI_SUM, 0, comm);
    if (my_pid == 0) printf("Glb #I %d\n", tot_num_I);
// #endif

    delete [] loc_map_f2c; delete [] glb_map_f2c;
    // MPI_Barrier(comm);
    // MPI_Abort(comm, -909);
    return C_rows_partits;
}
template int* coupled_pairwise(const par_CSRMatrix_I32F64F64_DOF1 *, par_CSRMatrix_I32F64F64_DOF1 * &, par_CSRMatrix_I32F64F64_DOF1 * &, const double );
template int* coupled_pairwise(const par_CSRMatrix_I32F64F64_DOF2 *, par_CSRMatrix_I32F64F64_DOF1 * &, par_CSRMatrix_I32F64F64_DOF1 * &, const double );
template int* coupled_pairwise(const par_CSRMatrix_I32F64F64_DOF3 *, par_CSRMatrix_I32F64F64_DOF1 * &, par_CSRMatrix_I32F64F64_DOF1 * &, const double );
template int* coupled_pairwise(const par_CSRMatrix_I32F64F64_DOF4 *, par_CSRMatrix_I32F64F64_DOF1 * &, par_CSRMatrix_I32F64F64_DOF1 * &, const double );

template int* coupled_pairwise(const par_CSRMatrix_I32F32F32_DOF1 *, par_CSRMatrix_I32F32F32_DOF1 * &, par_CSRMatrix_I32F32F32_DOF1 * &, const float  );
template int* coupled_pairwise(const par_CSRMatrix_I32F32F32_DOF2 *, par_CSRMatrix_I32F32F32_DOF1 * &, par_CSRMatrix_I32F32F32_DOF1 * &, const float  );
template int* coupled_pairwise(const par_CSRMatrix_I32F32F32_DOF3 *, par_CSRMatrix_I32F32F32_DOF1 * &, par_CSRMatrix_I32F32F32_DOF1 * &, const float  );
template int* coupled_pairwise(const par_CSRMatrix_I32F32F32_DOF4 *, par_CSRMatrix_I32F32F32_DOF1 * &, par_CSRMatrix_I32F32F32_DOF1 * &, const float  );

